%% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function out = Model_Schafer(CalibrationFactors,cfg,exvars,depvars,indepvars,r,t)
%% Destinguish between calibration or validation and specific type of run
global globals % Used to keep track of earlier iteration and for realtime plots
switch cfg.OutputSwitch
    case {'Normal','Normalt0'}
        
    case 'Calibration'
        indepvars.PreferenceFactor(r,t,squeeze(indepvars.UseModes(r,1,:))' & cfg.CalibrateOverModes)...
            = CalibrationFactors;
        
    case 'OptimizeTimeWeight'
        depvars.TimeWeight(r,t)          = CalibrationFactors; 
        
    case 'CalibrationPreferenceFactor0TimeWeight0'
        indepvars.PreferenceFactor(r,t,squeeze(indepvars.UseModes(r,1,:))' & cfg.CalibrateOverModes)...
            = CalibrationFactors(1:end-1);
        depvars.TimeWeight(r,t)          = CalibrationFactors(end); 
    otherwise
        disp('not a valid OutputSwitch')
end

indepvars.PreferenceFactor(r,t,~cfg.CalibrateOverModes)     = cfg.DefaultPreferenceFactor;
depvars.TravelPrice(r,t,:)                                  = exvars.TravelPrice(r,t,:); % Imported instead of vehicle module\
modes                                                        = indepvars.UseModes(r,1,:);

%% Model body, switch between different models
    switch cfg.Model{1}
        case 'Girod et. al, TIMER'
            %% Girod et. al TIMER model
            if ~(isfield(cfg,'SkipTimeWeight') && cfg.SkipTimeWeight == 1)
            depvars.TimeWeight(r,t) =...
                depvars.TimeWeight(r,t-1) .* sum(depvars.TimeUse(r,t-1,modes )...
                .* depvars.TravelDemand(r,t-1,modes ),3) ...
                ./ indepvars.TTB(t);
            end
    
        depvars.RelativeModeCost(r,t,modes ) =...
            depvars.TravelPrice(r,t,modes ) ./ mean(depvars.TravelPrice(r,t,modes )) .* indepvars.PreferenceFactor(r,t,modes )...
            + depvars.TimeUse(r,t,modes )  ./ mean(depvars.TimeUse(r,t,modes)).* depvars.TimeWeight(r,t);
                        
        depvars.ModeSplitOptimum(r,t,modes ) =...
                exp(-indepvars.Lambda .* depvars.RelativeModeCost(r,t,modes )) ...
                / sum(exp(-indepvars.Lambda .* depvars.RelativeModeCost(r,t,modes )));
       
        case 'Kyle and KIM, GCAM'
        %% Kyle and Kim, GCAM Model
        depvars.ModeSplitOptimum(r,t,modes ) =...
            indepvars.PreferenceFactor(r,t,modes) .* depvars.TravelPrice(r,t,modes).^indepvars.Lambda ./...
            sum(indepvars.PreferenceFactor(r,t,modes) .* depvars.TravelPrice(r,t,modes).^indepvars.Lambda);
    end
        
    if ~(isfield(cfg,'SkipModeSplit') && cfg.SkipModeSplit == 1)
    depvars.ModeSplit(r,t,modes )= ...
        depvars.ModeSplit(r,t-1,modes ) + indepvars.Inertia(modes ) .*( depvars.ModeSplitOptimum(r,t,modes ) - depvars.ModeSplit(r,t-1,modes ) ) ;
    elseif cfg.FixStartingModeSplit == 1
        depvars.ModeSplit(r,t,modes ) = exvars.TravelDemand(r,t,modes ) / sum(exvars.TravelDemand(r,t,modes ));
    else
        depvars.ModeSplit(r,t,modes ) = depvars.ModeSplitOptimum(r,t,modes );
    end
    
    depvars.TMB(r,t)             = 0.035 + 0.085 * sum(depvars.ModeSplit(r,t,5:end),3); 

    depvars.TravelDemand(r,t,modes) = ...
       depvars.ModeSplit(r,t,modes ) .* depvars.TMB(r,t) .* exvars.Income (r,t) ./ sum(depvars.ModeSplit(r,t,modes ) .* depvars.TravelPrice(r,t,modes),3);
   % Formulation above is equal to saying ModeSplit.* TMB .* Income ./mean( ModeSplit TravelPrice );
   
%% Write output
    depvars.ModeWeight(r,t,modes)= sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,:)),3);% Ensures not all modes are equal, but that car is not the only important mode
    
    switch cfg.OptimizationMethod{1}
        case {'Optimize over TimeWeight TTB'}

        case {'Iterative improvement of TimeWeight, optimization over PreferenceFactor','Iterative improvement of TimeWeight'}

    end
        
    switch cfg.OutputSwitch % Switches between possible output for different switches
        case {'Normal','Normalt0'}
            depvars.Resid(r,t)          = sum((depvars.TravelDemand(r,t,modes) - ...
            exvars.TravelDemand(r,t,modes)).^2./...
            exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes)).^2;
            out = depvars;
        case {'OptimizeTimeWeight'}
            depvars.Resid(r,t)          = (sum(depvars.TravelDemand(r,t,modes).*depvars.TimeUse(r,t,modes)) - ...
            indepvars.TTB(r,t)).^2./indepvars.TTB(r,t)^2;           
            out = depvars.Resid(r,t);
        case {'Calibration','CalibrationPreferenceFactor0TimeWeight0'}
            depvars.Resid(r,t)          = sum((depvars.TravelDemand(r,t,modes) - ...
            exvars.TravelDemand(r,t,modes)).^2./...
            exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes)).^2;
            out = depvars.Resid(r,t);
            globals.Error(1,end+1)              = out;
            globals.PreferenceFactor(:,end+1)   = squeeze(indepvars.PreferenceFactor(r,t,:)) ; 
            globals.TimeWeight(:,end+1)         = depvars.TimeWeight(r,t,1);        

        otherwise
            disp('not a valid OutputSwitch');
    end
        %% Visualizes minimization process
        plot_travel('LiveMinimizeVisuals','save',cfg,exvars,depvars,indepvars,output,r,t);
end